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ABSTRACT 



^ ' Adaptive finite elements are the method of choice for accurate simulations of optical components. However as 
Q shown recently by Bienstman et al. many finite element mode solvers fail to compute the propagation constant's 
imaginary part of a leaky waveguide with sufficient accuracy. In this paper we show that with a special goal 



^ , oriented error estimator for capturing radiation losses this problem is overcome. 
£^ ' Keywords: finite elements, goal oriented adaptive methods, PML, leaky modes 

Ph! 1. INTRODUCTION 



Radiation losses play an important role in the design of open resonator cavities and optical waveguides. Typically 
the field is strongly confined within such a resonant optical device. A prominent example is a leaky waveguide 
QQ . as discussed by Bienstman et al.""^ Figure 1 shows the corresponding waveguide geometry. The waveguide core 
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Figure 1. Leaky waveguide as discussed by Bienstman et al.^ The substrate-spacer-air layer stack is considered as 
infinite. Transparent boundary conditions are needed along the whole boundary. 



is inoTintcd on a spacer between the substrate and air. The electromagnetic field of the waveguide mode is 
strongly confined in a region around the waveguide core, Figure 2. Within the air and the spacer the field decays 
exponentially whereas the field propagates in the substrate. This wave propagation causes the radiation losses. 
The field amplitude is several orders of magnitude smaller within the substrate block than in the core. Therefore 
the imaginary part of the propagation constant is much smaller than the real part. For a precise computation of 
the propagation constant's imaginary part a sufficient accurate simulation of the electromagnetic field within the 
substrate is needed. Furthermore the transparent bomidary condition needs to be adapted. Standard adaptive 
finite element methods adapt the mesh such that the numerical solutions on successive meshes converge with 
optimal order in the energy norm. That is, the mesh is adapted such that the error 

Sh = [ (E - E^)e(E - E^) + (H - H^)/.(H - H^) 
Jn 

is minimized within the computational domain. Here (E, H) is the electromagnetic field and (E/j,Hfi) its finite 

element approximation. The real part of the propagation constant shows an optimal convergence as well. Since 
the imaginary part is small it may happen that the imaginary part is not improved at all when the mesh is 
refined. This is not a failure of the adaptive mesh refinement strategy because the error estimator was not 
designed for capturing the radiation losses. Hence one should apply a goal-oriented error estimator as proposed 
by Rannacher et al.'^ for various applications in engineering. This concept is based on a goal functional j(E) 
(quantity of interest). A special mesh refinement strategy is used to minimize the error 

S = |j(E)-j(E,.)|. 

In the Section 6 we will show that a special error estimator for capturing radiation losses allows for a robust 
and precise simulation of the imaginary part of the propagation constant. Before we will survey the finite 
element method and the goal oriented error estimation concept in the sections 2- 5. Functional expressions for 
the radiation losses which serve as the goal functional are derived in the Sections 3 and 5. We also address the 
construction of goal oriented error estimators for light scattering simulations. 

2. MAXWELL'S EQUATION ON UNBOUNDED DOMAINS AND FINITE 

ELEMENT DISCRETIZATION 

We start from the time-harmonic Maxwell's equations for the electric field 

Vx /x"^V X E-w^eE = -iu>3 
eV- E = 0. 



Here the electric field E and the material tensors are considered to be defined on the infinite space R"^. In 
this paper we restrict ourselves to waveguide geometries which exhibit a spatial invariance in the z— direction. 
Furthermore, the electromagnetic field as well the source current J should depend harmonically on z, 

B{x,y,z) = E{x,y)e''''^ 
Jix,y,z) = J{x,y)e''''^ 

with a propagation constant kz- In the following we drop the wiggles to keep the notation simple. With the 
definition Vfe^ = {d^, dy, kz)'^ the time-harmonic Maxwell's equations now read as 

Vfc^ X /x"^Vfe, X E(a;,y) -a;^eE(a;,y) = -iuj3{x,y) (la) 
eVfc, • E(a;,y) = 0. (lb) 

For the finite element method to be applicable it is necessary to bring Maxwell's equations into a variational 
form. To this end one needs to define a state space of admissible fields E. One encounters the difficulty that 
the fields of interest E have typically infinite field energies on the whole space R^. For example, a propagating 



Figure 2. Electric field's absolute value of the fundamental propagation mode for a waveguide as depicted in Figure 1 
with H — 220nm, D = Ifim and W = 0.5fim. The field is strongly confined in a region near the waveguide core. However, 
in this pseudocolored plot one perceives a refiection-free (no amplitude modulations) wave propagation in the substrate. 




Figure 3. Triangulation of the computational domain fl with PML. In the exterior domain we introduce a special 77) — 
coordinate system such that no material jump occurs in ^—direction. The boundaries of the quadrilaterals corresponds to 
isolines of the {^,ri) coordinate system. The exterior domain is truncated at ^ = p whichs yields a bounded domain Sip. 



mode of a loaky waveguide is exponentially increasing in certain directions towards infinity. Furthermore, above 
Maxwell's equations lack a boundary condition at infinity so far. Hence one could aim to define a state space which 
incorporates the outward radiation condition and yields integrable expressions in the variational formulation of 
Maxwell's equations. This is indeed done in the Pole condition approach, see.^' However, in this paper we use 
the perfectly matched layer method (PML). The PML method was originally invented by Berenger^ and analyzed 
and improved by several other authors. The essential idea is to replace the field E outside a computational 
domain by a complex continuation chosen such that one observes an exponential decay in the outward 
direction. The complex continued field E b satisfies Maxwell's equations with modified material tensors /U and e 
outside the computational domain. Due to its exponential decay it is allowed to truncate Eb to a finite domain 
rip. Here, p is a discretization parameter which corresponds to the thickness of the artificial sponge layer Qp \ D,, 
c.f. Figure 3. It can be shown that the truncation error disappears exponentially fast with increasing p for certain 
types of geometries.''"^ To deal with inhomogeneous domains it is necessary to define a special coordinate system 
(t], ^) in the exterior domain to guarantee the analyticity of E in the outward direction ^. Especially no material 
jump is allowed in the ^— direction, see paper. ^° To adjust the PML thickness we use an automatic strategy as 
proposed in our paper. The variational form is now given as 

/ Vfe, X * • pt"^Vfe, X E - • IE = -iw / ¥-J e Hk,{^p,cml) 

Jn„ Jn 



or in short notation 



a(*,E) = /(*) V*Gi/fc,(f^p,curl) 



with accordingly defined bilinear form a{-, ■) : Hk^ {^p^ cwcVjxHk^ (f^p, curl) C and fimctional / : Hk^ (ilp, curl) ^ 
C. Here Hk^{flp, curl) is the space of fields E(a;, y) G C"^ with finite electric and magnetic field energy within fip. 
The material tensors are modified in the exterior domain to incorporate the PML. The precise definition of the 
modified material tensors Jj, and e^in the exterior domain ilp \ are given in our paper. We again drop the 
wiggles for the sake of a simpler notation. In the following we will frequently make use of the formula 

/ Vfe, X E-pJ-^Vfe, X E- w^E-^E + iw*- J = -29 (fc^) / [E x /i"^Vfe, x E] ■ 
Jn Jq 

+ / [E X /z-^Vfc^ X E] -fi^y. (2) 
Jon 

Here = (0,0,1)''" and fixy is the normal on dil perpendicular to the z— direction. Formula (2) may be 
derived from Maxwell's equations by a multiplication of equation (la) with E{x,y) and integration over the 
three dimensional domain O x [0, zg] followed by partial integration based on Green's formula and afterwards 
differentiation with respect to zg. 

The variational form of Maxwell's equation (2) is discretized with high-order vector elements 'S'/i G C 
iJfc^ (S7p, curl). Here Vh is a finite dimensional subspace of iJ/c^ (fip, curl) spanned by the finite element ansatz 
functions •, *iv. This gives the discretized variational problem for E/j e V/j, 

a(*^,E;,) = fi^h) y^h^Vh, 

which can be casted into a linear algebraic system 

(Ao + fc.Ai + fcf As - w^B) g = f (3) 

Here e is the coefficient vector e = (ei, . . . , ejv) of the finite element solution. 
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3. RESONANCE MODES, WAVEGUIDE MODES AND SCATTERING SOLUTIONS 

Maxwell's equations give rise to three typical simulation scenarios which we shortly review below. In each case 
one computes certain quantities of interest from the simulated field E, e.g. the electric field energy Ea(E) within 
fl, the power flux Pn(E) through the cross section computational domain fl and the in plane power flux Pan(E) 
across the boundary of the cross section computational domain: 



Resonance Modes 



In this setting one seeks eigenpairs (E, cj) for given propagation constant and J = 0. Numerically this 
corresponds to an eigenvalue problem for (e, w^) in the matrix equation (3). We want to relate the imaginary 
part of u) to the functionals above. Assuming transparent materials in the computational domain fl and Q{kz) = 
one derives from formula (2) that 

Hence up to normalization the imaginary part of the resonance mode lo is proportional to the in-plane losses of 
the resonator. Therefore construction of a goal orientated error estimator for the computation of the qualtity 
factor of a cavity should start from (4). 

Waveguide Modes 

In this setting one seeks pairs (E, kz) for a prescribed angular frequency w G R such that equation (2) is satisfied 
with J = 0. Now as follows from equation (2) the imaginary part of kz is given by 

So the imaginary part of kz is exactly given as the ratio of the in plane power flux across the boundary of the 
computational domain and the power flux within the computational domain in waveguide direction. 

Scattering Solutions 

In this setting u) as well as kz are prescribed and either the source term J is non-zero or an incident fleld is 

present. It is described in-'^-'^ how to incorporate an incident field into the variational form of Maxwell's equation. 
Here we restrict ourselves to the special case that the current source consists of a "line" dipole source 

3{x, y, z) = S{x -xo,y- yo)e''''^3o 

with (a;o,yo) in ^- In this case the regularity of the solution is poor and the finite element method will suffer 
from a slow convergence. To cure that we use the subtraction approach, where the total field is split into an 
analytically given singular field Eg and a correction field Ec with smoothness properties fitting well in the finite 
element concept. We assume that the line source is not placed on a material jump, so that = ^i{xo,yo) and 
eo — e{xo,yo) are well defined. The solution Eg to Maxwell's equation 

Vfc^ X fi^'^Vk.x Es(x,y) - cj^eoEs(a;, y) = -iujd{x - xo,y - yo)Jo 

with constant material tensors /xq and eo is known analytically.^^ Now the correction field satisfies Maxwell's 
equation with the right hand side term 



fc = 



-Vfc, X (n ^-Mo^)Vfc, X Es(a;,y) + (e - eo)Es. 



Since fjr^ — fj,Q^ and e — eo are zero at (a;o,?yo) it can be proven that the right hand side fc is sufficiently 
smooth for an optimal convergence rate of the finite element method. The right hand side fc might be non-zero 
everywhere especially in the exterior domain. Due to the analyticity of fc the additional source terms will cause 
no problem in the PML. Since many finite element implementations do not allow for a source term in the PML 
one alternatively may compute the total field in the exterior domain. In this case additional boundary terms on 
dO, stemming from the singular field are needed. This is similar to the treatment of coupling with an incoming 
field and is naturally implemented within the finite element method. 



4. GOAL ORIENTED ERROR ESTIMATOR: GENERAL CONCEPT 

We follow the general concept for the construction of goal oriented error estimators as proposed by Rannacher.^ 
The quantity we aim to compute is a functional j(E) of the solution field E. So the finite clement mesh should be 
adapted such that the error _7(E) — j(E/,) is efficiently reduced. The idea by Rannacher is to embed this situation 
into the framework of optimal control theory. To this end, we define the trivial constrained optimization problem 

j(E)-j(E/.)= ^ min {jX*) - j(E/.) : a (*,*) = /(*) V* G iffe, (l^p, curl)} . (6) 

The Minimum of (6) corresponds to a stationary point (E, E*) of the Lagrangian 

£(E, E*) = i(E) - + f{E*) -a{E*, E) 

where E* denotes the 'dual' variable (Lagrangian multiplier). Hence, we seek the solution (E, E*) to the Euler- 
Lagrange system 

a(*,E) = /(*) G HkA^p,cml) (7a) 

a(E*,*) = /(E;*) V* G if^, (fip, curl). (7b) 

The first equation (7a) is the variational form of original Maxwell's equations. In the dual equation (7b) the 
functional j(E) - the quantity of interest - enters the stage in form of its linearization j'(E). A finite element 
discretization of the Euler-Lagrange system yields a supplemental discrete problem 

Introducing the primal and dual residuums, 



p(E^;.) = /(•)-a(-,E) 
piK;.) = j'{Ef,--)-a{El,-) 



one obtains the error representation 



i(E) - i(E^) = lp(Eft; E* - Wh) + ^p*(E^; E - Vh) + R (8) 

with remainder term R of third order in E — E/^ and E* — E^, see.^ 
Eigenvalue problems of the form 

a(*,E) - (t6(*,E) = 

can also be treated within this approach. To do this one extends the state space by the vector space of complex 
numbers C, V ^ V x C. Since a solution of an eigenvalue problem is only fixed uniquely up to normalization 
one includes a normalization, e.g. 6(E, E) = 1 into the variational form (9), 

a(*,E)-c76(*,E)+x(6(E,E)-l) = 0. V(*,x) e Ffc,(J7^,curl) x C. 

which gives the non-linear variational problem 



S((*,x),(E,a)) = 0. 



The corresponding dual problem is derived in."'^'' For a chosen goal functional j(E), it is necessary to check 
the unique solvability of the dual problem and its discretization by finite elements. For a wide range of goal 
functionals an error bound of the form 

b-(E-Eh|<C^|H|T||p1|T (9) 

T 

can be derived from the error representation (8), see.^ Here ||p||t denotes the norm of the residuum restricted 
to the triangle T. As usual for an adaptive strategy we refine triangles with large local contributions ||p||t||p*||t 
to the sum, see.^^ 



5. GOAL ORIENTED ERROR ESTIMATOR: ADAPTIVE COMPUTATION OF 

RADIATION LOSSES 

In this section we derive a goal oriented error estimator for the precise computation of radiation losses in optical 
devices. Here we restrict ourselves to waveguide mode problems. So the non-linear functional (5) is a natural 
candidate for the goal fmic;tional. But since the nominator on the left hand side of equation (5) involves derivatives 
of the solution field E this functional is not defined on the whole space Hk^{^lp, cml). A regularization of the 
functional is required. This is done in our case by an appropriate extension of the boundary integral expression 
into a volume integral around the boundary dQ, 

^ ^ (/r X M"^Vfc. X ^] • + X M"^Vfc. x ^] • ny) , . 

In our implementation we simplify this expression by replacing Vfe^ x * by Vfe^ x Eh- The so modified functional 
is unchanged when applied to the solution field E. 



6. A WAVEGUIDE EXAMPLE 

We apply our finite element solver to the challenging example of a leaky waveguide as depicted in Figure 1. This 
example was discussed on the OWTNM 2006 by Bienstman et al.^ and it was shown that an accurate computation 
of the losses is challenging. We computed the propagating mode with three mesh refinement strategies; uniform 
mesh refinement, adaptive mesh refinement with energy norm error estimator and adaptive mesh refinement with 
the new loss functional based error estimator. Third order finite elements were used. The computed eigenmodes 
are given in the Tables 1-3. In the first column the refinement steps are given, in the second the number of used 
unknowns and in the third column the computed effective refractive index. The reference value for the effective 
index of the waveguide is taken from,^ 

rieff = kjko = 2.4123720 + 2.91348 • 10"^^. (11) 

To start the eigenmode solver an initial guess for the propagation constant is needed. Such a guess can be 
obtained by placing a dipole source in the center of the waveguide for various kz, see Figure 4. Time-Harmonic 
Maxwell's equations are solved and the right hand side of Formula (5) is computed. In the resonant regime we 
observe a dramatic reduction of in-plane losses. The cusp of the curve in Figure 4 is an excellent initial guess for 
the eigenmode solver. 

In all three cases one observes convergence to the reference value. The benefits of an adaptive mesh refinement 
strategy are significant. As expected with the loss based mesh adaption (Table 3) the imaginary part of the 
propagation constant converges faster than with the energy based mesh adaption, especially for lower accuracy 
demands. As can be seen from Table 2 the imaginary part is not improved in the first three refinement steps 
when the energy based error estimator is used. However, for high accuracy demands they asymptotically show 
nearby the same behavior. For the loss based adaptive refinement method it takes less than 10 seconds to reach 
a relative accuracy of four digits in the imaginary part on a standard PC. 



Step 


N" DOF 


5R(nefr) 







3075 


2.40362444e+00 


1.76283e-08 


1 


6900 


2.41244790e+00 


2.95880e-08 


2 


16530 


2.41258532e+00 


2.90962e-08 


3 


43710 


2.41245837e+00 


2.91179e-08 


4 


129750 


2.41240459e+00 


2.91284e-08 


5 


428550 


2.41238429e+00 


2.91324e-08 


6 


1533030 


2.41237665e+00 


2.91339e-08 



Table 1. Computed effective refractive index: uniform mesh refinement. 



Step 


N" DOF 


5R(ncff) 









3075 


2.40362444e+00 


1 


.762836-08 


1 


3531 


2.41244233e+00 


1 


.570486-08 


2 


4581 


2.41257875C+00 


1 


.53788C-08 


3 


7080 


2.41246067e+00 


1 


.594416-08 


4 


12780 


2.41240706e+00 


2 


.960016-08 


5 


21918 


2.41238542C+00 


2, 


.960670-08 


6 


40875 


2.41237710e+00 


2 


.913686-08 


7 


73977 


2.41237393e+00 


2 


.913706-08 


8 


140631 


2.41237272C+00 


2 


.91348C-08 


9 


253395 


2.41237226e+00 


2 


.913486-08 


10 


476511 


2.41237209e+00 


2, 


.913486-08 



Table 2. Computed effective refractive index: energy based mesh refinement. In the first refinement steps the imaginary 
part is not improved. 



Step 


N° DOF 


5R(noff) 







3075 


2.40362444C+00 


1.76283(5-08 


1 


4698 


2.40636085C+00 


3.10891C-08 


2 


8895 


2.41259520C+00 


2.90980C-08 


3 


14718 


2.41250075e-h00 


2.911016-08 


4 


23019 


2.412425716-^00 


2.912406-08 


5 


39903 


2.41239249e-h00 


2.913076-08 


6 


62856 


2.41237984e-h00 


2.913336-08 


7 


112818 


2.41237497e-F00 


2.913436-08 


8 


203637 


2.41237312e-h00 


2.913466-08 


9 


385551 


2.41237241e-h00 


2.913486-08 



Table 3. Computed effective refractive index: loss based mesh refinement. Imaginary part is improved from the beginning 
of the mesh refinement process. 
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Figure 4. Ratio of the in plane power flux across the boundary of the computational domain and the power flux within the 
computational domain in waveguide direction. Time harmonic Maxwell's equations are solved for a vacuum wavelength 
Ao = 1.55/[im and various kz- As excitation a dipole is placed in the center of the waveguide core. One observe a dramatic 
reduction of in-plane losses when kz reaches the real part of the propagation mode. 
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